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Abstract 


In this paper, we introduce two new schemes based on the global Hessen- 
berg processes for computing approximate solutions to low-rank Sylvester 
tensor equations. We first construct bases for the matrix and extended 
matrix Krylov subspaces by applying the global and extended global Hes- 
senberg processes. Then the initial problem is projected into the matrix 
or extended matrix Krylov subspaces with small dimensions. The reduced 
Sylvester tensor equation obtained by the projection methods can be solved 
by using a recursive blocked algorithm. Furthermore, we present the upper 
bounds for the residual tensors without requiring the computation of the 
approximate solutions in any iteration. Finally, we illustrate the perfor- 
mance of the proposed methods with some numerical examples. 
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1 Introduction 


Let [y,1o,...,Iy © N. The multidimensional array ¥ = (Xisi5.-in) (1 < 
i; < 1j,j =1,...,N) is called an N-mode tensors with I, J2---In entries. 
There has been increasing research on tensors in recent years. For instance, 
Chang, Pearson, and Zhang [8] generalized the Perron—Frobenius theorem for 
nonnegative matrices to the nonnegative tensors. Eigenvalues, eigenvectors, 
symmetric hyperdeterminants were defined by Qi [31] for the real super- 
symmetric tensors, and their properties were described. In [30], the restart 
techniques are described for the tensor infinite Arnoldi method. 


In this work, we introduce two new projection methods for solving the 
low-rank Sylvester tensor equation 


My AY $i eg AS orc ba Ky AY = B, (1) 


where the matrices AM € RI», n = 1,2,...,N, and right-hand side 
tensor B € R2*2* IN are given, and X € R2*2**IN is an unknown 
tensor. The Sylvester tensor equation (1) has a unique solution if and only if 
Ap tAg+:::+Aw £0, for all A; € o0(AM), i=1,2,...,N, where o(A™) is 
the spectral of matrix AM [9]. In this study, it is assumed that the Sylvester 
tensor equation has a unique solution. The Sylvester tensor equations are 
one of the famous problems arising from the discretization of a linear partial 
differential equation in high dimensions by the use of finite elements, finite 
differences, and spectral methods [27, 28, 37]. The Sylvester matrix equation 

AYX + XAQ?T — B, 
is a special case of the Sylvester tensor equation (1), where X is a 2-mode 
tensor. Many iteration methods for computing approximate solutions for 
the Sylvester tensor equations (1) have been introduced in recent years. For 
example, Chen and Lu [9] proposed the GMRES method based on tensor 
form (GMRES-BTF) to solve the Sylvester tensor equation. Also, to speed 
up the convergence of the GMRES-BTF method, they proposed precondi- 
tioned GMRES-BTF. Beik, Saberi Movahed, and Ahmadi-Asl [4] presented 
some iterative methods based on the tensor format to solve the Sylvester 
tensor equations (1). In [33, 34], Saberi-Movahed et al. introduced the ten- 
sor format of restarted Simpler GMRES, (SGMRES-BTF(m)), to solve the 
Sylvester tensor equation and described an accelerating method in accor- 
dance with a modification of the generalized conjugate residual with inner 
orthogonalization (GCRO) method based on the tensor format. Bi-conjugate 
gradient (BiCG) and bi-conjugate residual (BiCR) methods as well as their 
preconditioned versions based on the tensor format, have been presented in 
[39]. The tensor form of the global least squares method is proposed in [24]. 
Huang, Xie, and Ma [22] proposed the tensor form of the GMRES method 
for solving a class of tensor equations via the Einstein product. Furthermore, 
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for the case in which the coefficient tensor is symmetric, they proposed the 
MINRES and SYMMLQ methods based on the tensor format. Dehdezi and 
Karimi [15] extended the conjugate gradient squared and the conjugate resid- 
ual squared methods to solve the generalized coupled Sylvester tensor equa- 
tions. In [16], the authors proposed a gradient based iterative method version 
for solving the tensor equations and presented a new preconditioner to accel- 
erate the convergence rate of the proposed iterative methods. A projection 
method has been introduced in [3] to find approximations of linear systems 
in low-rank tensor format. Kressner and Tobler [25] proposed the Krylov 
subspace for the case in which the right-hand side tensor has a low-rank. 
Recently, Bentbib, El-Halouy, and Sadek [5] introduced a new projection 
method to compute approximate solutions for the low-rank Sylvester tensor 
equations. The extended Krylov-like methods were proposed in [6] to find the 
solutions for the low-rank Sylvester and Stein tensor equations. The block 
and extended block Hessenberg algorithms for solving the Sylvester tensor 
equation with low-rank right-hand side (1) were presented in [12]. Hessen- 
berg based methods are among the popular methods in terms of the Krylov 
subspace methods, with less need for arithmetic operations and less storage 
space compared to the Arnoldi-based methods. The Hessenberg process con- 
structs nonorthogonal bases for the associated Krylov subspace. ‘The schemes 
based on the Hessenberg process have recently received great attention; see, 
for instance, [32, 35, 19, 17, 21, 12]. This motivated us to introduce two new 
projection schemes, employing the global Hessenberg process on the matrix 
Krylov subspaces. The main idea of this scheme is to project the problem 
onto a matrix or an extended matrix Krylov subspace. Then the reduced 
problem can be solved by using the recursive blocked algorithm [11]. Com- 
plexity consideration is given to show that the global and extended global 
Hessenberg processes are less expensive than the global and extended global 
Arnoldi ones. 


We use the following notations. For the matrices X and Y in R”*”, 
we consider the following inner product (X,Y) = tr(X7Y), where tr(-) 
denotes the trace. The associated norm is the Frobenius norm denoted by 
||E\|p. The notation XL-Y means that (X,Y)r = 0. The n x n identity 
matrix is denoted by I‘. Moreover, ey denotes the 7th canonical vector 
of R®, and 0imxn denotes the m x n zero matrix. 


The remainder of this paper is organized as follows. In section 2, we 
review some basic notations and definitions. In section 3, the global Hessen- 
berg process with maximum strategy and an approach for solving (1) with 
a right-hand side tensor of a specific rank is described. The extended global 
Hessenberg approach is presented in section 4. The complexity of the new 
methods is considered in section 5. Some numerical examples for evaluating 
the performance of our approaches are given in section 6. Finally, section 7 
gives a brief conclusion. 
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2 Preliminaries 


In this part, the notations and basic definitions of tensors are presented. 
Throughout this paper, we denote tensors by Euler script letters. Matrices 
and vectors are denoted by capital and lowercase letters, respectively. Also, 
the Kronecker product of matrices A and B is denoted by A ®@ B and the 
Kronecker product of tensors A and B, is denoted by A@B. Norm of an Nth 
order tensor ¥ € R!*!2*-"*4 ig denoted by ||4||~ and is defined as follows: 


qh Ig In 


Alle = V(X) =| D2 RG iy 


44=1lig=1 in=l 


Definition 1 ((13]). Denote the N-mode (matrix) product of a tensor V € 
RM Xf2x-xIn and a matrix U € RY*!" by X x, U. It is of dimension I, x 
Ig X +++ X In-1 X J X Ingi X +++ X In and defined as 


In 
(x Xn Octet = S Miz ig- in Ujin 
in=l 


Proposition 1 ((13]). Let A € R"*”*~*4» be an Nth order tensor, let 
BeR?*!m™,C € R*¥*!n, and let W € R&*!». Then 


A Xm BxnC=AXnC Xm B, 
A xn W XnC=AxX,CW. 


Definition 2 ([14]). Assume that ¥ € R“*/2*"*/» is an Nth order tensor 
and that {U} is a set of matrices U, € R’"*!"(n = 1,...,N). Then their 
product in all possible modes (n = 1,2,...,.N) is of size I, x In x +--+ x In 
and defined as follows: 


Xx {U} =X x1 U1 X2U2--- xn UN, 
and 
Xx {U}? =X x, U7? x2? > xy UE. 


Definition 3 ({13]). . The outer product of two tensors A € R'1*/2%- x1 
and B E€ Rix J2X- xX IN is denoted by Ao B E€ R11 Xf Int XS x Ja Xo XT 
with entries 


Cay tnedidw = Aa cin Bay sin : 


If vj, v2,...,un are N vectors of sizes I;,i = 1,...,N, then their outer prod- 
uct is an Nth order tensor of size , x Iz x --- x Iy and is given by 
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UL O° OUN’,..in = Ui(i1) ++ UN (in). 


Definition 4 ([{13]). An Nth order tensor ¥ € RD*!2*"*4n js called a rank 
one tensor if it can be written as the outer product of N vectors a; € R! (i = 
1,...,.N) as follows: 

X =a, 0d20-:-Oan. 


If a tensor can be written as a sum of R rank one tensors, then it is called a 
rank R tensor. 


Definition 5 ([26]). The Kronecker product of two tensor A = a10a20:--oan 
and B = bi 0 bg 0-:-0 by is defined as 


A® B= (a1 ® b1) 0-:- 0 (an @ by). 


Proposition 2 ([5]). Assume that A € R2*”2*~*!~ and B € R1X2* xiv 
are Nth order tensors, that A € R*»*!, and that B € R/™*/". Then 


(A @B) xn (A® B) =(A Xp A) @ (BX, B). 
Proposition 3 ((5]). The product of a rank one tensor A = a, 0a20---oay € 
R1Xf2x-xIn and a set of matrices U;, € R’”*!”,(n = 1,...,.N) is defined 
as follows: 


Ax {U} = Uja1 0--+-0Unan. (2) 


Definition 6 ((13]). The CP decomposition of an Nth order tensor A € 
Roxhx--xIn ig written as follows: 
R 
A= Soa) 0 a) o++-0 aN), 


r=1 


where R € N and a e€ R%, (i = 1,...,N). Assume that al, = 


1,...,N), are normalized. Then the CP decomposition is given by 
R 
A= yo rAaY oa o---oaW), 
r=1 
where A, € R. 


Definition 7 (Left inverse[35]). Consider Z, € R"** as a matrix partitioned 
as follows: 


Zk 


where Zi, is a k x k matrix. If the matrix Z,, is nonsingular, then a left 
inverse of Z;, is defined as follow 
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L —1 
ZK = [Zig 1 Onx(n—2)] - 


Definition 8 ((7]). Let A = [Ai, A2,...,A,] and B = [Bi, Bo,...,Bi] be 
matrices of dimension n x ps and n x Is, respectively, where A; and B; (i = 
1,...,p}j =1,...,0) are n x s matrices. Then the o-product of matrices A 
and B denoted by A’ © B is the p x | matrix defined by: 


(A? o B)ij = (Ai, Bj) r. 
Some properties that are verified by the ®- and products are as follows: 
1. (DA)? oR = AT e(D7 8). 
2. AT o(B(L@I)) = (AT oB)L. 


In what follows, we assume that the right-hand side B in (1) is of rank R. 
As known [13], by using the CP decomposition, 6 can be written as 


R 
B= Sb o--- bf), (3) 


where BO = ee an sh | € R&X® j = 1,...,N, are the factor ma- 


trices. By simple calculations, we can rewrite the relation (3) as 
B=Trx, BY... xy BO), (4) 


in which Zp denotes the identity tensor of Nth order of size Rx Rx---x R 
with ones along the super-diagonal. 


3 Global Hessenberg process with maximum strategy 


The global Hessenberg process with maximum strategy was first presented 
in [17] by Heyouni to build a basis of the matrix Krylov subspace 


m—-1 
Km(A,V) = {So nay where 7; € R for f=0tenmaah 
i=0 


where A € R”*” and V € R”*®*. The global Hessenberg process with maxi- 
mum strategy can be summarized in Algorithm 1 [17]. 

By employing Algorithm 1 with m = m; and s = R for the pair 
A®, BO), we obtain Vinig1 = [VO?,..., Vey] € ROM D2 with VO € 


R°*®. for k =1,...,m,;+1, and the upper Hessenberg matrix H,,, = (ni) € 
Rim+1)xmi which satisfy 


AOV rag = Vint (Am; ® i), (5) 


IJNAO, Vol. 12, No. 3 (Special Issue), 2022, pp 658-679 


664 Cheraghzadeh, Toutounian and Khoshsiar Ghaziani 


Algorithm 1 The Global Hessenberg process with Maximum Strategy 


1. Input: Nonsingular matrix A, initial block V, and an integer m. 
2. Determine ip and jo such that |Vig,jo| = max {|Vi,j|}1 225° 5 B=Vioior 
Vi =V/B; ly =%0; c1 = jo. 


3. For k=1,2,...,m 

4 U = AV,. 

5. For j=1,2,...,k 

6. hy = U1; ,0;3 U=U—hy pV;. 

7 End For. 

8 Determine ig and jo such that |Ui,,j;5| = max {|Uj,;|)}} S752 
hrti,e = Vig,jos Vert = U/hr4ijes let+1 = 40; Ceti = Jo- 


9. End For. 


= Vin; (Am, ® 7) 1 hae Ve ee ® 7), (6) 


a 


where H,,, denotes the matrix obtained from Fe by deleting its last row. 
As [5], we consider an approximate solution of (1) as 


Xm = (Vn @ Tr) x {Vin}, (7) 


where {V,,,} denotes a set of matrices {Vm,,VWm.,---,Wmy} and Vp, is an 
my, X-+++xX my tensor satisfying the low-dimensional Sylvester tensor equation 


N 
ye Ym xi Am, = BEm, (8) 
i=1 

where § = [Jj 6; and Em = (e{"") 0--- oe"). 


Proposition 4. Let R,, be the residual tensor corresponding to the approx- 
imate solution ¥,, of (1). Then 


N 
Rm >= ES Am +1,mi (Vin Xi elms)” 


i=l 


) @ Ie X1 Ving Xi Ve XN Vins 


(9) 


where Y,, is the solution to (8). 


Proof. The proof is similar to that of Proposition 6 in [12]. 


Theorem 1. Let 4, be an approximate solution of (1). Then the corre- 
sponding residual R,,, satisfies 


N 
m 
[Rll $ y/(@nR = (m= 1))2)% Do | Pacts |? [Yon Xi eR, II?, (10) 
i=1 


where m= max ™,. 
1<i<N 
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Proof. The proof is similar to that of Theorem 7 in [12]. 


Furthermore, from the fact that 
[VnjI? <nmjR, 6=1,...,N, 


we have 


? [Yn Xa ein, II? (11) 


N 
|Rml| < \/ (nmR)N ys | Pmi+dmi 
i=1 


The upper bounds (10) and (11) are pessimistic. We propose the following 
approximation, which is derived heuristically, 


N 
|Rin|| © Em = V/(nmR) a | Pengtlymi 


i=1 


2 [Yin Xi Chr, 


2. (12) 


Similar to Algorithm 2 in [5], the global Hessenberg process with the max- 
imum strategy for the Sylvester tensor equation (1) can be summarized in 
Algorithm 2. 


Algorithm 2 


1. Input: Coefficient matrices A“) ,i = 1,...N, and the right-hand side in low-rank 
representation, 
B= (Bo, Be), = eae) . 


2. Output: An approximate solution Ym for equation (1). 

3. Choose a tolerance « > 0, integer parameters ki, i=1,...,N. Set kj =0,m; = ky. 

4. Fori=1:N 

6. Construct the basis [Vita ake Vin! | and the matrix Hm, by 
Algorithm 1. . 

t: End For 

8. End For 

9. Solve the low-dimensional equation x. Ym Xi Hm; = BEm by the recursive 


blocked algorithms presented in [11]. 
10. Compute the estimated residual norm of Rm, 


i.e, Bm = X/(ranR)y/N, | Rmetims |? [Ym Xi eB, |. 
11. If Em >, set kj = kj 4 ki, m, = ki 4 k, for i=1,...,N, and go to step 4. 
12. Compute the approximate solution by %m = (Ym @Z™) x1 Vy + XN Vins 
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4 The extended global Hessenberg process 


We first recall the extended matrix Krylov subspace. Let A € R”*” and 
V € R”**. The extended global Hesssenberg process corresponding to the 
pair (A,V) is defined as follows [17]: 


Ke,(A,V) = span(V, A~!V, AV,...,A"™-2V, A-™V), 
=Ky,(A,V)+ | onn  aee A-tYy), 


The algorithm proceeds by running one step of the Global Hessenberg process 


with A and one step with A~!, while maintaining orthogonalization among 


Es 
all generated vectors and the n x s matrices Y; = eeten whose entries 


are zero except (Yj)1;,c; = 1. The first two block vectors Vi and Ve are 
obtained as follows: 


ye =V/ru, (13) 


<j< 
where 131 = Vj, ,c, and |Vi, ¢,| = max{|Vjj|}iten and 


Ve?) = W/ro, 25 (14) 


where W = AV —ri2V, r112= (ATV) i 5 72,2 = Wis,ca and |Wi.c9| = 
<j<s 
max{|Wi,j|}1fen- 
Let V; = [V,?, V,)] be the ith n x 2s block vector of Vm = [Vi,-.-,Vml 
and let 
hai-1,25-1 hai-1,25 
Ay = , 
haiaj-1—hai,ag 
be the 2 x 2 block matrix (i,7) of the upper block Hessenberg matrix H,, € 
R2™+1)x2m_ Then we compute the two block vectors Ve and Ve) by the 
relation 


k 


se ve (Heti. @ T°) = [Av 1- SJ] V,V, il (Hyp ® Pr), 
j=l 

(15) 

where the entries of coefficients matrices Hy+1,, and H;.,, for i = 1,...,k, 


will be determined such that the relations 
1 1 
Ve leVisecceYor and (VA) Vee = 1, 


and e (2) 
Veg len, wee) YOR and (1555 beeieenes = 
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hold for k = 1,...,m. The determination of indices lo441,co~n41 and 
lop+2,Cok+2 is similar to that of indices 1;,c, and lg,c2, respectively. The 
main steps of the extended global Hessenberg process algorithm to generate 
V., and H,,, may be summarized as follows. 


Algorithm 3 The Extended Global Hessenberg process with Maximum 
Strategy 


1. Input: Nonsingular matrix A, initial block V, and an integer m. 

2. Determine ip and jo such that |Vig,j9| = max {|Vj,j|}} S758 ; 1,1 = Vio,jo3 
Ve? =V/ria; L =%t0; c1 = jo; - 

3. W=A-!V; T1,2 = Wii cy: 


1 1<j< 
4,.W=W-maVi, — |Wigiol = max{|Wi,s I} SIS%; 72,2 = Wig,jo3 


v2) = W/r2,2; lg=%to, c2=Jo- 


5. For k=1,2,...,m 
6. Wwe= AV), 
7. For i=1,...,k 
1 
: hoi-1,2k-1 = Wig; 1,ca3-15 W=W- ee ), 
2 
hoi2k—1 = Wig;,co.. W=W- nA ), 
9. End For. . 
10. Determine i9 and jo such that |Wi,,j9| = max {|Wi,j|}7 525° ; 
hop+1,2k-1 = Wig, joi Vey = W/hor+1,2k-1; loe+1 = 403 Con+1 = jo- 
11. W= A-1V 2), 
12. For i=1,...,k. 
ie hai—1,2k = Wig;_1,c21-11 W = W—- hoi—1,26V5; 
2 
h2iak = Wig;,co,33 W = W — haianV5 2 
14 End For. 
1 
15 hor+1j2k = Wrong 1 cons? W = W — hon st,26V)- 
16. Determine i9 and jo such that |Wi,,j| = max {|Wi,j|}H S753 ; 


2 : . 
har+2,2k = Wig, 503 Any = W/hor+2,2k3 lokt+2 = 40; Cor+2 = Jo- 
17. End For. 


Suppose that the matrix P,,, is defined by [Y1, Y2,..., Yam]. Then 
Fo Ve = Le, 


where L,,, € R?”*?”™ is a unit lower triangular matrix. So, we have L;,)(P7,,° 
Vm;) = 1@™), As in [1], we consider V4, = (Pm(Lz7 @ I&)))T = (LZ) @ 
I‘s))P? | as a left inverse for the product, which verifies the relation V4, o 
Vm = 1?™s), Using this matrix, we can state the following proposition. 


Proposition 5. Let T,, = Veg © (AV,,), and suppose that m steps of 
Algorithm 3 have been carried out. Then 


AV = Vent Tm OT), (16) 
= Vin (Tn ed I()) + Visa eae & T°), (17) 
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where T;,; is the 2 x 2 block (i,j) of Tm and EX = [02x2¢m—1), 1], and Tm 
is obtained by removing the two last rows of T,». 


Proof. The proof is similar to the case for the classical Arnoldi process in 
[20]. 


As [86], in the following proposition, we derive some recursive relations, 
which can be used to significantly reduce the computational cost of the basic 
algorithm. 


Proposition 6. Let Ti; = [t.1,...,¢:2m] and Hy», = [h.1,...,h:,2m] be two 
2(m +1) x 2m block upper Hessenberg matrices, let (+) = (¢;;) = Hyp 
and let 71,1, 71,2,72,2 be as defined in Algorithm 3. Then for the odd columns, 
we have 

t. 25-1 = h: 25-1, j = 1, ey MM, 


and for the even columns, we have 


1 m 
(k = 1) t.2 = aries vt) 
2°93 
tan (ser — [ baa he. 
(2m—2) x2 


— Ti et:1), 


2)\ 1 (2 
1 = By, 
(l<k<m) ton = tox + t:on-10™, 


2(m+1 TA: k+1 
fanaa = (cpr — [Tehran |, 


pe) = (te. 
Proof. Starting from (15), we have 
AV? = Vieur (Hev1nel? ® 1°) + Ve (Hey, @ 1) 
= Vier (Hee, & Fe"), 


Pre-multiplying the above relation by V4, 41, we get 


Vint1 ° AV = Vint1 ° VieiHaese @ 1) 
=> (VP ag ° View1) Hee”, 


[(2k+2) ro (9k) 
= He 
ite ns eres poked 
shoal? 
Omani 


Hence, 
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t.or—1 = hs 2e-1, k=1,...,m 
From the lines 2 and 3 of Algorithm 3, we have 
r2V,?) = nr AVY) - ravi. 
Pre-multiplying this relation by A, we get 
12.9 AV,”) = ravi? — 2 AV,”. 


Pre-multiplying the above relation by V%, 41, we have 
1 
(VE4.0 AV?) = = (ria(Vinga 9 VO yr Vig OAV). 


Consequently, 


In addition, from (15), one gets 
VP = AVer1(Hasanes” @ 1) + AVe (Heese @ I). 
This relation implies that 
VE 41 0 AVisi (Hii nes” @ 1) 
= VE o Vi) — ” 410 (AV, (Hye @ 1?) 
a = et) (Ving © AVE Hes,” 


) 
= er(m+1) _ | Teh1:2k,2k I: 
ae O(2m—2k) x 2k 


On the other hand, for the left-hand side of this relation, we deduce 
VE 1° AVeg1 (Hes ney? @ I) 


" 1) (2) hors1261 
= Vat [AV AV 41 bee 


= hope1,2kW 41 © AVE + hop+2,2kW 41 © AVE, 


= hogs+i2rt:or41 + hor+e2,2Kt:2h42- 


Hence 


1 


cote Tyht:2k,2k 
are hok+2,2k 


(—hor41,2Kt:,2h41 + enn: f 
(2m—2k) x 2k 
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By using the inverse of the 2 x 2 upper triangular matrix H,+41,, and defining 
A= (ERY) 1 pd) this relation can be written as follows: 


7 Tht: 
ae = trope pet) fe (e2( +1) | k/1:2k,2k uae 


O(2m—2k) x 2k) 


which completes the proof. 


4.1 Extended global Hessenberg process for low-rank 
Sylvester tensor equation 


In this subsection, we consider the extended global Hessenberg process de- 
rived in the previous subsection for the pair (AM, B®), i = 1,...,N. By 
applying Algorithm 3 with s = R to the pair (A,B), i=1,...,N, the 
block matrices V,,, = Ves = 1,...,N, are obtained and the 


following relation holds, fori =1,...,N, 
AV mn; = Meet (Ties ® pe) 


SVT Ol EV TN eg BEI), 518) 


+1,mi 


where ED = [02x2,...,02x2,1)] € R?X?™, and Tra, = (Th?) € R20 +D) x2m: 
is the restriction of AM to the extended global Krylov subspace K¢, (AM, BM). 
Using Line 1 of Algorithm 3, we have 


BO =r VP), for += 1,2,...,N. 


As in the case of the global Hessenberg process, for the low-rank Sylvester 
tensor equation (1), we seek an approximate solution of the form 


where {V,,} denotes a set of matrices Vn, € R"*2?™, 7 = 1,...,N, and 

Vm € R?2™*"*2™N satisfies the low-dimensional Sylvester tensor equation 
N 
> Yin x4 Tm; = BmEm, (20) 
i=l 

where Bm = TI, r{) and Em = (e&?™? o-+-0e?™"), In this case, the 


residual corresponding to ¥,, can be written as 
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N 


Rei re SET BO TG Viren REV er Ri Vin, 
= (21) 
We can easily obtain 
N 
Real] < of (n= 2m + LYm)N | Son Xi TO). m BZ, (22) 


w=1 


and 


[Rll < 1 (QnmR)% 3 [Yn Xs Tita a EF, I (23) 


w=1 


where m = max Mi. Finally, the following estimate is derived heuristically: 
StS 


[Rmll ~ Em = \/(2nmR) mp? |Yn x fu, Mt, pee o 


(24) 


For the extended global Hessenberg process, the main part of Algorithm 
2 remains the same except that the lines 6, 9, and 10 must be changed as 
follows: 
6. Construct the basis [Vist is dg Ve. 4K and the matrix T,,, by Algo- 
rithm 3 and the formulas of Proposition 6. 


9. Solve the low-dimensional equation pal Ym Xi Tm; = BmEm by the 
recursive blocked algorithms presented in [11]. 


10. pits the estimated residual norm of R,,, that is, 


= V/GnmR)/ ON 1 |lYm x acne ee 


2. 


5 Complexity consideration 


In this section, we present the required number of operations to solve the 
low-rank Sylvester tensor equation (1) for 1) = Ig =--- = In. Let Nnz de- 
note the number of nonzero elements of matrix A, and suppose that the LU 
decomposition of A is available for computing the block matrix W = A7!V. 
We compare the required operations for the extended global Hessenberg pro- 
cess and the extended global Arnoldi process [18]. oe 3 requires 
(2n?s +4ns) operations for computing the block matrices Ve? and V2. In 
addition, the iteration k of this algorithm involves 
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e Vas which requires 2s Nnz + ns(4k + 1) — 4k? operations, 


° Vee which requires 2n?s + ns(4k + 3) — (2k + 1)? operations. 


Note that the global Arnoldi process (Algorithm 2 in [18]) requires 2n?s+10ns 
operations for computing the global QR decomposition [V, A~!V], and the 
iteration k of this process involves 


eU= [av,™, A-1y,), which requires 2sNnz + 2n?s operations. 


e Hy; =VJ oU, U=U—-V,(Hij @IT°), i=1,2,...,k, which require 
16nsk operations. 


e the global decomposition of U, that is, U = Viri(Hpsi., ® 1), which 
requires 10ns operations. 


Therefore, for computing an approximation of the solution of Sylvester tensor 
equation (1), the total cost number of operations required to perform m 
iterations of the extended global versions of Arnoldi and Hessenberg processes 
is approximately shown in Table 1. In addition, the total cost number of 
operations required to perform m iterations of the global Hessenberg process 
(Algorithm 1) and the modified global Arnoldi process (Algorithm 2.2 in 
[23]) is presented in this table. According to Table 1, when solving the low- 
rank Sylvester tensor equation (1), the global and extended global Hessenberg 
processes are less expensive than the global and extended global Arnoldi ones. 
On the other hand, these Hessenberg processes use the maximum strategy. 
Hence they involve some data movement. However, these processes need 
slightly less storage than the Arnoldi processes per iteration. 


Table 1: Operation count for the global and extended global versions of 
Hessenberg and Arnoldi processes. 


Process Number of operations 

Global Arnoldi N(QmRNnz + (m+ 1)(2m + 3)nR — (m(m + 1))/2) 

Global Hessenberg N(2mRNnz + (m+ 1)?nR — (m(m + 1)\(2m+ i ))/6) 

Extended Global Arnoldi N(2mRNnz + 2(m + 1)n?R + (m+ 1)(8m + 10)nR) 

Extended Global Hessenberg N(2mRNnz + 2(m+1)n?R+4(m + 1)?nR — m(8m? + 18m + 13)/3) 


6 Numerical experiments 


In this section, some test problems with N = 3 are used to examine the ro- 
bustness of two new presented methods for solving the low-rank Sylvester 
equation (1). All the numerical experiments were performed in double- 
precision floating-point arithmetic in MATLAB 2021a. The machine we have 
used is an Intel(R) Xeon(R) CPU E5-2680 v4@2.40 GHz, 128 GB of RAM, 
using the Tensor Toolbox [2]. We employ the recursive blocked algorithms 
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introduced in [11] to solve the low-dimensional Sylvester tensor equations (8) 
and (20). The step size parameter k’ associated with one cycle is equal to 3. 
The algorithms stopped whenever E,, < 10~’, where E,, is the estimate of 
||Fm||. We also compare the numerical behavior of the methods in terms of 
the number of cycles (Cycle), the norm of residual ||R ||, the norm of error 
||X** — X,,||, where ** is the exact solution, and the CPU time in seconds 
(CPU time) required only for constructing the Krylov subspace basis and the 
solution of reduced Sylvester tensor equation. Note that we use the proce- 
dure cp_als(B, R) from the toolbox [2] to compute the CP decomposition of 
the right-hand side B. In Table 2, we report ||B — B-p||, where the Bey is the 
CP decomposition corresponding to the right-hand side tensor B, using the 
procedure cp_als(B, R). The results of examples are reported in Table 2. For 
each example, the rank R and the dimension n are presented in this table. 
In Figure 1 , by plotting the norm of residual || R,,||~” versus the number of 
cycles, we display the convergence history of the global and extended global 
Arnoldi and Hessenberg algorithms for Examples 1-5. 


Example 1. In this example, as in [5], we consider the matrices AM,i = 
1, 2,3, corresponding to discretization of the operator 
Ou 


0 
Lu) = Au- fle war + blew) z+ ole) 


in the unit square [0,1] x [0,1] with Dirichlet homogeneous boundary condi- 
tions. The number of inner grid points in each direction is no for the operator 
L. The discretization of the operator DL yields matrices extracted from the 
Lyapack package [29], using the command fdm and denoted as 


A = fdm(no, fi(x,y), fo(a,y), g(a, y)), = 1, 2,3, 


with fi(,y) _ ee foley) = sin(x,y), g(x,y) = y? ~ o = no. The 
right-hand side tensor is chosen in such a way that the exact solution of the 
Sylvester tensor equation (1) has the form ¥* = x, 0 22043, with a; = 
rand(n, 1), for ¢ = 1, 2,3. 


Example 2. Assume that in the Sylvester tensor equation (1), the coefficient 
matrices are presented as [5] 


A = gallery('poisson’,no), i= 1,2,3, 


where n = n2. The right-hand side tensor is constructed such that the exact 
solution ¥ of the Sylvester tensor equation (1) is a tensor with entries equal 
to one. 


Example 3. Let A, i = 1,2,3, be defined as [10] 


A® = rand(n,n) + diag(ones(n,1) * alfa), 
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where alfa = 8 and the right-hand side tensor is constructed as in Example 
1. 


Example 4. Consider the Sylvester equation (1) with the coefficient matrices 
generated by [38] 


A® = diag(rand(n-1,1),-1) + diag(2 + diag(rand(n,n))), i = 1,2,3, 


and the right-hand side tensor is constructed as in Example 1. 


Example 5. The coefficient matrices AM, i = 1, 2,3, for the Sylvester tensor 
equation (1) are defined as 


ij : 1 
Al (1,9) = T+ (—J]’ 


and the right-hand side tensor is constructed as in Example 1. 


Table 2: Results of Examples 1-5. 


Example Algorithm |B -— Bep|l Rll |’* — Xn || Cycle CPU time 
Global Arnoldi 3.655e —08 8.549e-—08 9.908e — 11 30 2.879 
Example 1 Global Hessenberg 3.655e —08 2.90le—07 2.667e — 10 28 2.558 
n= 400 ,R=4 Extended Global Arnoldi 3.655e-—08 4.197e-—08 3.173e—11 7 0.261 
Extended Global Hessenberg 3.655e—08 1.41le—O07 1.162e—10 6 0.110 
Global Arnoldi -355e—08 1.406e—08 1.560e — 08 14 0.138 
Example 2 Global Hessenberg -355e —08 1.573e—08  1.735e — 08 14 0.229 
n = 400 , R=3 | Extended Global Arnoldi .395e—08 1.375e—08  1.608e — 08 5 0.079 
Extended Global Hessenberg 1.355e—08  4.528e€—08 2.652e — 08 4 0.058 
Global Arnoldi 532e—-05 1.53le—05  3.731le — 07 19 0.612 
Example 3 Global Hessenberg 532e—05 1.530e-—05  3.729e — 07 18 0.479 
n= 500 , R=3 | Extended Global Arnoldi .532e—-05 1.53le—05 3.73le — 07 9 0.429 
Extended Global Hessenberg 1.532e—05 1.53le—05  3.73le — 07 8 0.267 
Global Arnoldi -980e —O7 1.980e—O07 2.698 — 08 5 0.049 
Example 4 Global Hessenberg -980e —07 1.984e—-O07 2.704e — 08 4 0.046 
n=500 .R=3 Extended Global Arnoldi 980e—O07 1.980e—O7 2.698e — 08 3 0.082 
Extended Global Hessenberg 1.980e—07 1.980e—O07 2.698e — 08 3 0.077 
Global Arnoldi .038e —08 1.034e—08 2.567e — 09 12 0.120 
Example 5 Global Hessenberg -038e —08 1.16le—08 2.622e — 09 1 0.144 
n= 500 , R=3 | Extended Global Arnoldi 1.038e —08 1.042e—08 2.567e — 09 5 0.115 
Extended Global Hessenberg 1.038e—08  1.034e—08 2.566e — 09 5 0.112 
As can be seen from Table 2 and Figure 1, Global Arnoldi, Extended 
Global Arnoldi, and Global Hessenberg, Extended Global Hessenberg meth- 


ods are shown a similar behavior. In addition, for all examples, the number 
of cycles of Extended Global Hessenberg is less than or equal to that of the 
other methods. In Examples 1, 2, 3, and 5, the CPU time of Extended 
Global Hessenberg method is less than the others. The results of Example 4 
show that when the required number of cycles is small for Global Hessenberg 
method, this method outperforms the other methods in terms of CPU times. 
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Figure 1: Convergence history of the global and extended global Arnoldi and 
Hessenberg algorithms for Examples 1-5. 


7 Conclusion 


In this study, for computing the approximate solutions of the Sylvester tensor 
equation (1) with the low-rank right-hand side, two new projection methods 
based on the Hessenberg process were proposed. The theoretical results of 
these methods were presented and analyzed as well. The global and extended 
global Hessenberg algorithms were compared, in terms of CPU times, cycles, 
and the number of operations, with the global and extended global Arnoldi 
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algorithms, respectively. Numerical examples showed that the global and 
extended global Hessenberg algorithms are efficient and feasible for solving 
the low-rank Sylvester tensor equation (1). 
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